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Abstract 

Adaptive mesh refinement teciiniques are nowadays an establislied and 
powerful tool for the numerical discretization of PDE's. In recent years, 
wavelet bases have been proposed as an alternative to these techniques. The 
main motivation for the use of such bases in this context is their good per- 
formances in data compression and the approximation theoretic foundations 
which allow to analyze and optimize these performances. We shall discuss 
these theoretical foundations, as well as one of the approaches which has been 
followed in developing efficient adaptive wavelet solvers. We shall also discuss 
the similarities and differences between wavelet methods and adaptive mesh 
refinement. 
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1. Introduction 

Among those relevant phenomenons which are modelled by partial differential 
or integral equations, countless are the instances where the mathematical solutions 
exhibit singularities. Perhaps the most classical examples are elliptic equations 
on domains with re-entrant corners, or nonlinear hyperbolic systems of conserva- 
tion laws. While such singularities are sources of obvious theoretical difficulties — 
classical solutions should be abandonned to the profit of weak solutions — they are 
also an obstacle to the convergence of numerical approximation methods, in the 
sense that they deteriorate the rate of decay of the error with respect to the size of 
the discrete problem : achieving a prescribed accuracy will typically require finer 
resolution and therefore heavier computational cost and memory storage, in com- 
parison to the approximation of smooth solutions. Let us remark that singularities 
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often have a physical relevance : they represent the concentration of stress in elas- 
ticity, boundary layers in viscous fluid flows, shock waves in gas dynamics... It is 
therefore a legitimous requirement that they should be accurately resolved by the 
numerical method. 

In this context, the use of adaptive methods, appears as a natural solution 
to improve the approximation at a reasonable computational cost. Here, the word 
adaptivity has a twofold meaning : (i) the discretization is allowed to be refined only 
locally, in particular near the singularities of the solution, and (ii) the resolution 
algorithm uses information gained during a given stage of the computation in order 
to derive a new refined discretization for the next stage. The most typical example 
is adaptive m,esh refinement based on a-posteriori error estimates in the finite ele- 
ment context. While these methods have proved to be computationally successful, 
the theory describing their advantages over their non-adaptive counterpart is far 
from being complete. In particular, the rate of convergence of the adaptive algo- 
rithm, which describes the trade-off between the accuracy and complexity of the 
approximation, is not clearly understood. 

In recent years, wavelet bases have been proposed as an alternative to adaptive 
mesh refinement, motivated by their good performances in data (more specifically 
image) compression. In wavelet-based adaptive schemes, the set of basis functions 
which describe the approximate solution is updated at each stage of the compu- 
tation. Intuitively, the selection of the appropriate basis functions plays a similar 
role as the selection of the mesh points in adaptive finite element methods, and 
one could therefore expect similar performances from both approaches. On a more 
rigorous level, a specific feature of the wavelet approach is the emergence of a sound 
theoretical setting which allows to tackle fundational questions such as the rate of 
convergence of the adaptive method. 

The goal of this paper is to give some elements of comparison between adaptive 
wavelet and mesh refinement methods from this perspective. We shall first describe 
in §2 a general setting which leads us in §3 to a first comparison between wavelets 
and adaptive finite elements from the point of view of approximation theory. We 
discuss in §4 the relation between these results and adaptive algorithms for PDE's. 
After recalling in §5 the classical approach in the finite element context, we present 
in §6 an adaptive wavelet strategy which has been applied to various problems, and 
discuss its fundational specificities. Finally, we shall conclude in §7 by pointing out 
some intrinsic shortcoming of wavelet-based adaptive methods. 



2. A general framework 

Approximation theory is the branch of mathematics which studies the process 
of approximating general functions by simple functions such as polynomials, finite 
elements or Fourier series. It plays therefore a central role in the accuracy analysis of 
numerical methods. Numerous problems of approximation theory have in common 
the following general setting : we are given a family of subspaces {Sn)n>o of a 
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normed space X, and for / G X, we consider the best approximation error 

a^if) := inf (1) 

TypicaUy, N represents the number of parameters which are needed to describe an 
element in Sn, and in most cases of interest, ffNif) goes to zero as this number 
tends to infinity. If in addition (JnH) < CN^^ for some s > 0, we say that / is 
approximated at rate s. 

Given such a setting, the central problem of approximation theory is to char- 
acterize by some analytic (typically smoothness) condition those functions / which 
are approximated at some prescribed rate s > 0. Another important problem is 
how to design simple approximation procedures / /at € Eat which avoid solving 
the minimization problem while remaining near optimal in the sense that 

<Ca^(/), (2) 

for some constant C independent of N and /. 

As an example, consider approximation by finite element spaces Vu defined 
from regular conforming partitions Th of a domain Vl C into simplices with 
uniform mesh size h. The approximation theory for such spaces is quite classical, 
see e.g. and can be summarized in the following way. If denotes the 

classical Sobolev space, consisting of those functions / e such that -D"/ G V 
for I a I < t, we typically have 

/ e W'+^-P inf 11/ - g\\w^., < Ch' (3) 

provided that Vh is contained in W^'^ and that Vh has approximation order larger 
than t + r, i.e. contains all polynomials of degree strictly less than t + r. Such 
classical results also hold for fractional smoothness. We can express them in terms 
of the number of parameters, remarking that := dim(V;i) ~ h^"^, so that if we 
set X = W^'P and Sn ■= Vh with h := N'^^'^, we have obtained 

f eW'+''^P ^ aNif) <CN~''/'^. (4) 

We have thus identified an analytic condition which ensures the rate s = r/d. 
Note that this is not a characterization (we only have an implication) , yet a deeper 
analysis shows that an "if and only if" result holds if we slightly modify the notion 



of Sobolev smoothness (using Besov classes, see 13 ). In summary, the rate of 
approximation in W'P is governed by the approximation order of the Vh spaces, 
the dimension d and the level of smoothness of / measured in L^. Let us finally 
remark that near-optimal approximation procedures can be obtained if we can find 
a sequence of finite element projectors Pjv : X i— * Sn such that ||PAr||j(:-»x < K 
with K independent of A^ : in this case, we simply take /at = -P/v/ and remark that 
\\f-fN\\x < il+K)aNif). 

In the following we shall address the same questions in the cases of adaptive 
finite element and wavelet approximation. As we shall see, a specific feature to such 
cases is that the spaces Sat are not linear vector spaces. 
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3. Adaptive finite elements and wavelets 

In the adaptive finite element setting, the number of parameters N is pro- 
portional to the number of triangles, but for a given budget N the partition T 
and the finite element space Vr are allowed to be locally refined in a way which 
depends on the function / to be approximated. It is therefore natural to define the 
approximation spaces Sn as 

Sn ■= ^#{T)<nVt- (5) 

It should be well understood that the Sn are not linear vector spaces (the sum 
of two elements does not in general fall in Sn when their triangulation do not 
match) but any g e Sn is still described by 0{N) parameters, which encode both 
its triangulation T and its coordinates in Vr- The requirement of adaptivity has 
thus led us to the concept of nonlinear approximation. 

Wavelet bases offer another track toward nonlinear adaptive approximation. 
The simplest prototype of a wavelet basis is the Haar system. Let us describe this 
system in the case of expanding a function / defined on [0, 1] : the first component 
in this expansion is simply the average of /, i.e. the orthogonal projection (/, eo)eo 
onto the function cq = X[o,i]- The approximation is then refined into the average 
of / on the two half intervals of equal size. This refinement amounts in adding the 
orthogonal projection (/, ei)ei onto the function Ci = X[o,i/2] — X[i/2,i]- Iterating 
this refinement process, we see that the next components have the same form as ei 
up to a change of scale : at refinement level j, we are adding the projection onto 
the functions 

^jMx) = Vl'^i,[Vx -k), fc = 0, • • ■ , 2^' - 1, (6) 

where ip — ei. Since all these functions are orthogonal to the previous ones, letting 
j go to +CX), we obtain the expansion of / into an orthonormal basis of L'^{[0, 1]) 

/ = E /^V'A, (7) 

AGV 

with fx := {f,ip\). In the above notation A concatenates the scale and space 
parameters j and k, and V is the set of all indices (including also the first function 
eo). In order to keep track of the scale j corresponding to an index A = {j,k) 
we shall use the notation |A| — j. More general wavelet systems in one or several 
space dimension are built from similar nested approximation processes, involving 
e.g. spline functions or finite elements in place of piecewise constant functions (see 
or for a general presentation). 

This brief description suggests that a natural construction of adaptive wavelet 
approximations is obtained by using only a limited set of indices A as the scale |A| 
grows, which depends on the function to be approximated and typically corresponds 
to those wavelets whose supports are close to its singularities. It is therefore natural 
to define the approximation spaces Sn as the set of all TV terms combinations 



Sn :={5]dA^A ; #(A) <iV}. 



(8) 
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Again this is obviously not a linear space, since we allow to approximate a function 
by choosing the best N terms which differ from one function to another. Note that 
we still do have Sn + Sn = S2n- 

Both adaptive finite element and wavelet framework have obvious similiarities. 
However, the answer to the two basic questions raised in the previous section — what 
are the properties of / which govern the decay of (Tjv(/) and how to compute in a 
simple way a near optimal approximation of / in Sjv — is only fully understood in 
the wavelet framework. Concerning the first question, a striking result by DeVore 
and his collaborators is the following : with X :— best TV-term wavelet 

approximation satisfies 

/ e W'+''"''' => <7Nif ) < CN-'^'^ (9) 

with q and r connected by the relation 1/q = 1/p + r/d, assuming that the mul- 
tiresolution approximation spaces associated to the wavelet basis are in W'^'P and 
contain the polynomials of degree strictly less than t + r. 

Such an estimate should be compared with the linear estimate (Q) : the same 
convergence rate is governed by a much weaker smoothness assumption on / since 
q < p (as in the linear case, an "iff and only if" result can be obtained up to 
slight technical modifications in the statement of (|^)). This result gives a precise 
mathematical meaning to the spatial adaptation properties of best iV-term wavelet 
approximation: a function / having isolated discontinuity, has usually a smaller 
amount of smoothness s + t when measured in L'' than when measured in with 
1/q = 1/p+t/d, and therefore cnH) might decrease significantly faster than ejv(/)- 

The answer to the second question is given by a result due to Temlyakov : if 
/ = X^Aev ^aV'Aj and if we measure the approximation error m. X = W^'^, a near 
optimal strategy when 1 < p < cxo consists in the thresholding procedure which 
retains the N largest contributions II^a^'aIIx ■ if A^r is the corresponding set of 
indices, one can prove that there exists C > independent of TV and / such that 

II/- ^ d^yjx\\x<CaN{f). (10) 

AeAiv 

This fact is obvious when X — using the orthonormal basis property. It is a 
remarkable property of wavelet bases that it also holds for more general function 
spaces. In summary, thresholding plays for best iV-term wavelet approximation an 
analogous role as projection for linear finite element approximation. 

In the adaptive finite element framework, a similar theory is far from being 
complete. Partial answers to the basic questions are available if one chooses to 
consider adaptive partitions with shape constraints in terms of a uniform bound on 
the aspect ratio of the elements 

max(^[Diam(ii:)]''/vol(X)) < C 

Such a restriction means that the local refinement is isotropic, in a similar way to 
wavelets. In such a case, we therefore expect a rate of approximation similar to 
(H). Such a result is not available, yet the following can be proved jl^ for Lagrange 



(11) 
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finite elements of degree m : if for any given tolerance e > 0, one is able to build 
a partition T = T(e) of cardinality N = N{£) such that on each K ^ T the local 
error of approximation by polynomials satisfies 



with q and r connected by the relation l/q — l/p + r/d and assuming s + t < m. 
The effective construction of T{e) is not always feasible, in particular due to the 
conformity constraints on the partition which does not allow to connect very coarse 
and very fine elements without intermediate grading. However, this result shows 
that from an intuitive point of view, the adaptive finite element counterpart to 
wavelet thresholding amounts in equilibrating the local error over the partition. One 
can actually use these ideas in order to obtain the estimate for adaptive finite 
elements under the more restrictive assumption that l/q < 1/p + r/d. Let us finally 
mention that the approximation theory for adaptive finite elements without shape 
constraints is an open problem. 

4. Nonlinear approximation and PDE's 

Nonlinear approximation theory has opened new lines of research on the theory 
of PDE's and their numerical discretization. On the one hand, it is worth revisiting 
the regularity theory of certain PDE's for which the solutions develop singularities 
but might possess significantly higher smoothness in the scale of function spaces 
which govern the rate of nonlinear approximation in a given norm than in the scale 
which govern the rate of linear approximation in the same norm. Results of this 
type have been proved in particular for elliptic problems on nonsmooth domains 
||2C| and for scalar ID conservation laws These results show that if u is the 

solution of such equations, the rate of decay of (Jn{u) is significantly higher for best 
iV-term approximation than for the projection on uniform finite element spaces, 
therefore advocating for the use of adaptive discretizations of such PDE's. 

On the other hand these results also provide with an ideal benchmark for 
adaptive discretizations of the equation, since ajqiu) represents the best accuracy 
which can be achieved by N parameters. In the wavelet case these parameters 
are typically the N largest coefficients of the exact solution u. However, in the 
practice of solving a PDE, these coefficients are not known, and neither is the set A 
corresponding to the indices of the N largest contributions ||c?a'0a||- It is therefore 
needed to develop appropriate adaptive resolution strategies as a substitute to the 
thresholding procedure. Such strategies aim at detecting the indices of the largest 
coefficients of the solutions and to compute them accurately, in a similar way that 
adaptive mesh refinement strategies aim at contructing the optimal mesh for finite 
element approximation. In both contexts, we could hope for an algorithm which 
builds approximations mat G Sjv such that ||m — UAr||x is bounded up to a fixed 




(12) 



then we can build global approximants G Vt C Sn such that 



(13) 
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multiplicative constant by aN{u) for a given norm of interest, but this requirement is 
so far out of reach. A more reasonable goal is that the adaptive strategy exhibits the 
optimal rate of approximation : if aN{u) < CN~^ for some s > 0, then — uatUx < 
CN^'^ up to a change in the constant. Another requirement is that the adaptive 
algorithm should be scalable, i.e. the number of elementary operations in order 
to compute uj^ remains proportional to N. Let us finally remark that the norm 
I • \\x for which error estimates can be obtained is often dictated by the nature of 
the equation (for example X = in the case of a second order elliptic problem) 
and that additional difficulties can be expected if one searches for estimates in a 
different norm. 



5. The classical approach 

The classical approach to numerically solving linear and nonlinear partial dif- 
ferential or integral equations T{u) = by the finite element method is typically 
concerned with the following issues : 

(cl) Well-posedness of the equation, i.e. existence, uniqueness and stability of the 
solution. 

(c2) Discretization into a finite element problem J-T-{uq-) = by the Galerkin 
method with uq- e Vr, analysis of well-posedness and of the approximation 
error \\u — ut||x- 

(c3) Numerical resolution of the finite dimensional system. 

(c4) Mesh refinement based on a-posteriori error estimators in the case of adaptive 
finite element methods. 

Several difficulties are associated to each of these steps. First of all, note that 
the well-posedness of the finite element problem is in general not a consequence 
of the well-posedness of the continuous problem. Typical examples even in the 
linear case are saddle point problems. For such problems, it is well known that, 
for Galerkin discretizations to be stable, the finite element spaces for the different 
solution components have to satisfy certain compatibility conditions (LBB or Inf- 
Sup condition), which are also crucial in the derivation of optimal error estimates. 
Thus the discrete problem does not necessarily inherit the "nice properties" of the 
original infinite dimensional problem. Concerning the numerical resolution of the 
discrete system, a typical source of trouble is its possible ill- conditioning, which 
interferes with the typical need to resort on iterative solvers in high dimension. An 
additional difficulty occuring in the case of integral equations is the manipulation 
of matrices which are densely populated. 

Finally, let us elaborate more on the adaptivity step. Since more than two 
decades, the understanding and practical realization of adaptive refinement schemes 
in a finite element context has been documented in numerous publications [Q, H, ^ 
^ Key ingredients in most adaptive algorithms are a-posteriori error esti- 
mators which are typically derived from the current residual J-{ut) ■ in the case 
where the Frechet derivative DJ-{u) is an isomorphism between Banach function 
spaces X to Y , one can hope to estimate the error ||m — UT-jlx by the evaluation of 



614 



Albert Cohen 



||-^(wT)||y- The rule of thumb is then to decompose ||jr(u7-)||Y into computable lo- 
cal error indicators rjn which aim to describe as accurately as possible the local error 
on each element K ^T. In the case of elliptic problems, these indicators typically 
consist of local residuals and other quantities such as jumps of derivatives across 
the interface between adjacent elements. A typical refinement algorithm will subdi- 
vide those elements K for which the error indicator rjK is larger than a prescribed 
tolerance e resulting in a new mesh T. Note that this strategy is theoretically in 
accordance with our remarks in §3 on adaptive finite element approximation, since 
it tends to equilibrate the local error. Two other frequently used strategies consist 
in refining a fixed proportion of the elements corresponding to the largest rjK, or 
the smallest number of elements K for which the r]K contribute to the global error 
up to a fixed proportion. It is therefore hoped that the iteration of this process 
from an initial mesh Tq will produce optimal meshes ('7^i)n>o in the sense that the 
associated solutions u„ := U7-„ G Vt„ converge to u at the optimal rate : 

cJn{u) < ^\\u^ ujx < C[#(r„)]-^ (14) 

up to a change in the constant C. Unfortunately, severe obstructions appear when 
trying to prove ( p^ even in the simplest model situations. One of them is that tjk is 
in general not an estimate by above of the local error, reducing the chances to derive 
the optimal rate. For most adaptive refinement algorithms, the theoretical situation 
is actually even worse in the sense that it cannot even be proved that the refinement 
step actually results in a reduction of the error by a fixed amount and that u„ 



converges to m as n grows. Only recently |26, 30| have proof of convergence appeared 



for certain type of adaptive finite element methods, yet without convergence rate 
and therefore no guaranteed advantage over their non-adaptive counterparts. 



6. A new paradigm 

Wavelet methods vary from finite element method in that they can be viewed 
as solving systems that are finite sections of one fixed infinite dimensional system 
corresponding to the discretization of the equation in the full basis. This observa- 
tion has led to a new paradigm which has been explored in ]l5| for linear variational 
problems. It aims at closely intertwining the analysis — discretization — solution pro- 
cess. The basic steps there read as follows : 

(nl) Well-posedness of the variational problem. 

(n2) Discretization into an equivalent infinite dimensional problem which is well 
posed in i'^. 

(n3) Devise an iterative scheme for the ^2-pi'oblem that exhibits a fixed error re- 
duction per iteration step. 

(n4) Numerical realization of the iterative scheme by means of an adaptive appli- 
cation of the involved infinite dimensional operators within some dynamically 
updated accuracy tolerances. 

Thus the starting point (nl) is the same. The main difference is that one aims at 
staying as long as possible with the infinite dimensional problem. Only at the very 
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end, when it conies to applying the operators in the ideal iteration scheme (n4), 
one enters the finite dimensional realm. However, the finite number of degrees of 
freedom is determined at each stage by the adaptive application of the operator, so 
that at no stage any specific trial space is fixed. 

The simplest example is provided by the Poisson equation — Au = / on a 
domain Q with homogeneous boundary conditions, for which the variational formu- 
lation 'm X ~ Hq reads : find u ^ X such that 

a{u,v) ~ L{v), for all v ^ X, (15) 

with a{u,v) := J^VuVv and L{v) := J^f'^- The well-posedness for a data / G 
X' = H^^ is ensured by the Lax-Milgram lemma. In the analysis of the wavelet 
discretization of this problem, we shall invoke the fact that wavelet bases provide 
norm equivalence for Sobolev spaces in terms of weighted l"^ norms of the coefficients 
: if u, = u\4>\, one has 

IM!1.-^-Eii-a^^iIh-'-E2'^'"I"^i'- (16) 

A A 

We refer to and ||2^ for the general mechanism allowing to derive these equiva- 
lences, in particular for Sobolev spaces on domains with boundary conditions such 
as Hq. Therefore, if we renormalize our system in such a way that HV'aUx = 1, we 
obtain the norm equivalence 

hlll-rf, (17) 

where U :— (uA)Aev and || • || denotes the £^ norm. By duality, one also easily 
obtains 

ll/llx' - \\Fr, (18) 
with F := ((/, ip\))\ev- The equivalent l'^ system is thus given by 

AU = F, (19) 

where A(A, /x) = a{ip\,'ipfj) is a symmetric positive definite matrix which is contin- 
uous and coercive in In this case, a converging infinite dimensional algorithm 
can simply be obtained by the Richardson iteration 

f/" [/"-I + t[F - ylC/"-i) (20) 

with < r < 2[Aniax(^)]~^ and — 0, which guarantees the reduction rate 
\\U - C/"|| < p\\U - ;7"-i|| with p = max{l - TAmin(A), rAmax(A) - 1}. Note that 
renormalizing the wavelet system plays the role of a multiscale preconditioning, 
similar to multigrid yet operated at the infinite dimensional level. 

At this stage, one enters finite dimensional adaptive computation by modifying 
the Richardson iteration up to a prescribed tolerance according to 

jjn jjn-i ^ T-(cOARSE(i^, e) - APPROX(AC/"-\ e)) (21) 

where ||i^ - COARSE(F, e) || < e and - APPROX(A[/"-i, £)|| < e, and the 
J7" are now finite dimensional vector supported by adaptive sets of indices A„. The 
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procedure COARSE, which simply corresponds to thresholding the data vector 
at a level corresponding to accuracy e, can be practically achieved without the 
full knowledge of F by using some a-priori bounds on the size of the coefficients 
(/i'0a), exploiting the local smoothness of / and the oscillation properties of the 
wavelets. The procedure APPROX deserves more attention. In order to limitate 
the spreading effect of the matrix A, one invokes its compressibility properties, 
namely the possibily to truncate it into a matrix An with N non-zero entries per 
rows and columns in such a way that 

\\A-AN\\f2^e2<CN~\ (22) 

The rate of compressibility s depends on the available a-priori estimates on the 
off-diagonal entries A(A,/i) :— Vtp\\/ipf_i which are consequences of the smooth- 
ness and vanishing moment properties of the wavelet system, see [Q. Once these 
properties are established, a first possibility is thus to choose 

APPROX(A[/"-\e) = AjvC/"-i (23) 

with N large enough so that accuracy e is ensured. Clearly the modified iteration 

( pl| ) satisfies ||C/ — C/"|| < p\\U — U"~^\\+2t£, and therefore ensures a fixed reduction 

rate until the error is of the order -r^e, or until the residual F — AU" is of order 

i-p ' 

^i-p^^ £. A natural idea is therefore to update dynamically the tolerance e, which is 
first set to 1 and divided by 2 each time the approximate residual COARSE(i^, e) — 
APPROX ( A?7"~ e) is below [ ^[i^^^ -l-3]£ (which is ensured to happen after a fixed 
number of steps). 

We therefore obtain a converging adaptive strategy, so far without information 
about the convergence rate. It turns out that the optimal convergence rate can also 
be proved, with a more careful tuning of the adaptive algorithm. Two additional 
ingredients are involved in this tuning. 

Firstly, the adaptive matrix vector multiplication APPROX has to be de- 
signed in a more elaborate way than (^3|) which could have the effect of inflating 
too much the sets A„. Instead, one defines for a finite length vector V 

3 

APPROX(At/, e)=J2 ^2^-' 1^2' - V21-1] (24) 

where Vat denotes the restriction of V to its iV largest components (with the notation 
^1/2 = 0)1 E^nd j is the smallest positive integer such that the residual X];=o ~ 
^23-! II II — V2i-i II -|- \\A\\ ||T^— V2J II is less than e. In this procedure, the spreading of 
the operator is more important on the largest coefficients which are less in number, 
resulting in a significant gain in the complexity of the outcome. 

Secondly, additional coarsening steps are needed in order to further limitate 
the spreading of the sets A„ and preserve the optimal rate of convergence. More 
precisely, the procedure COARSE is applied to C/" with a tolerance proportional 
to £, for those n such that e will be updated at the next iteration. 
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With such additional ingredients, it was proved in that the error has the 
optimal rate of decay in the sense that 

ctn{u) < CN-' ^ \\u - UnWx lit/ - C/"|l < C[#(A„)]-^ (25) 

and that moreover, the computational cost of producing u„ remains proportional to 
#(A„). It is interesting to note that this strategy extends to non-elliptic problems 
such as saddle-points problems, without the need for compatibility conditions, since 
one inherits the well-posedness of the continuous problem which allows to obtain a 
converging infinite dimensional iteration, such the Uzawa algorithm or a gradient 
descent applied to the least-square system (see also The extension to 

nonlinear variational problems, based on infinite dimensional relaxation or Newton 
iterations, has also been considered in It requires a specific procedure for 

the application of the nonlinear operator in the wavelet coefficients domain which 
generalizes (p^). It should also be mentioned that matrix compressibility also applies 
in the case of integral operators which have quasi-sparse wavelet discretizations. 
Therefore several of the obstructions from the classical approach — conditioning, 
compatibility, dense matrices — have disappeared in the wavelet approach. 

Let us finally mention that the coarsening steps are not really needed in the 
practical implementations of the adaptive wavelet method (for those problems which 
have been considered so far) which still does exhibit optimal convergence rate. 
However, we do not know how to prove ( |25| ) without these coarsening steps. There 
seems to be a similar situation in the finite element context : it has recently been 
proved in that ( p^ ) can be achieved by an adaptive mesh refinement algorithm 
which incorporates coarsening steps, while these steps are not needed in practice. 



7. Conclusions and shortcomings 

There exist other approaches for the development of efficient wavelet-based 
adaptive schemes. In particular, an substantial research activity has recently been 
devoted to multiresolution adaptive processing techniques, following the line of idea 
introduced in p8| , pof . In this approach, one starts from a classical and reliable 
scheme on a uniform grid (finite element, finite difference or finite volume) and 
applies a discrete multiresolution decomposition to the numerical data in order to 
compress the computational time and memory space while preserving the accuracy 
of the initial scheme. Here the adaptive sets A„ are therefore limited within the 
resolution level of the uniform grid where the classical scheme operates. This ap- 
proach seems more appropriate for hyperbolic initial value problems [ p^ p^ , in 
which a straightforward wavelet discretization might fail to converge properly. It 
should again be compared to its adaptive mesh refinement counterpart such as in 




Let us conclude by saying that despite its theoretical success, in the sense of 
achieving for certain classes of problems the optimal convergence rate with respect 
to the number of degrees of freedom, the wavelet-based approach to adaptive nu- 
merical simulation suffers from three major curses. 



618 



Albert Cohen 



The curse of geometry : while the construction of wavelet bases on a rectangu- 
lar domains is fairly simple — one can use tensor product techniques and inherit the 
simplicity of the univariate construction — it is by far less trivial for domains with 
complicated geometries. Several approaches have been proposed to deal with this 
situation, in particular domain decomposition into rectangular patches or hierarchi- 
cal finite element spaces, see [ p^ [2l} |, and concrete implementations are nowaday 
available, but they result in an unavoidable loss of structural simplicity in compar- 
ison to the basic Haar system of §3. 

The curse of data structure : encoding and manipulating the adaptive wavelet 
approximations C/" to the solution means that we both store the coefficients and 
the indices of the adaptive set A„ which should be dynamically updated. The same 
goes for the indices of the matrix A which are used in the matrix-vector algorithm 
( p^ ) at each step of the algorithm. This dynamical adaptation, which requires 
appropriate data structure, results in major overheads in the computational cost 
which are observed in practice : the numerical results in Q] reveal that while the 
wavelet adaptive algorithm indeed exhibits the optimal rate of convergence and 
slightly outperforms adaptive finite element algorithms from this perspective, the 
latter remains significantly more efficient from the point of view of computational 
time. 

The curse of anisotropy : adaptive wavelet approximation has roughly speaking 
the same properties as isotropic refinement. However, many instances of singu- 
larities such as boundary layers and shock waves, have anisotropic features which 
suggests that the refinement should be more pronounced in one particular direction. 
From a theoretical point of view, the following example illustrate the weakness of 
wavelet bases in this situation : if / = Xn with Q C M*^ a smooth domain, then 
the rate of best A^-term approximation in X — is limited to r = l/{2d — 2) and 
therefore deteriorates as the dimension grows. Wavelet bases should therefore be 
reconsidered if one wants to obtain better rates which take some advantage of the 
geometric smoothness of the curves of dicsontinuities. On the adaptive finite ele- 
ment side, anisotropic refinement has been considered and practically implemented, 
yet without a clean theory available for the design of an optimal mesh. 

The significance of wavelets in numerical analysis remains therefore tied to 
these curses and future breakthrough are to be expected once simple and appropriate 
solutions are proposed in order to deal with them. 



References 

[1] Babushka, I. and W. Reinhbolt (1978) A-posteriori analysis for adaptive finite 
element computations, SIAM J. Numer. Anal. 15, 736-754. 

[2] Babuska, I. and A. Miller (1987), A feedback finite element method with a- 
posteriori error estimation: Part I. The finite element method and some basic 
properties of the a-posteriori error estimator, Comput. Methods Appl. Mech. 



Adaptive Methods for PDE's Wavelets or Mesh Refinement? 



619 



Engrg. 61, 1-40. 

[3] Bank, R.E. and A. Wciscr (1985), Som,e a posteriori error estimates for elliptic 

partial differential equations, Math. Comp., 44, 283-301. 
[4] Barinka, A., T. Barsch, P. Charton, A. Cohen, S. Dahlke, W. Dahmen, K. 

Urban (1999), Adaptive wavelet schemes for elliptic problems — Implementation 
and numerical experiments, IGPM Report #173 RWTH Aachen, SIAM J. Sci. 
Comp. 23, 910-939. 

[5] Bcrgcr, M. and P. CoUcla (1989) Local adaptive mesh refinement for shock 

hydrodynamics, J. Comp. Phys. 82, 64-84. 
[6] Berger, M. and J. Ohger (1984) Adaptive mesh refinement for hyperbolic partial 

differential equations, J. Comp. Phys. 53, 482-512. 
[7] Bertoluzza, S. (1995) A posteriori error estimates for wavelet Galerkin methods, 

Appl. Math. Lett. 8, 1-6. 
[8] Bertoluzza, S. (1997) An adaptive collocation method based on interpolating 

wavelets, in Multiscale Wavelet Methods for PDEs, W. Dahmen, A. J. Kurdila, 

P. Oswald (eds.), Academic Press, 109-135. 
[9] Bihari, B. and A. Harten (1997) Multiresolution schemes for the numerical 

solution of 2-D conservation laws, SIAM .J. Sci. Comput. 18, 315-354. 
[10] Binev, P., W. Dahmen and R. DeVorc (2002), Adaptive finite element methods 

with convergence rates, preprint IGPM-RWTH Aachenm, to appear in Nu- 

mcrischc Mathcmatik. 
[11] Canuto, C. and I. Cravcro (1997), Wavelet-based adaptive m,ethods for 

advection- diffusion problems. Math. Mod. Mcth. Appl. Sci. 7, 265-289. 
[12] Ciarlet, P.G. (1991), Basic error estimates for the finite element method, Hand- 
book of Numerical Analysis, vol II, P. Ciarlet et J.-L. Lions eds., Elsevier, 

Amsterdam. 

[13] Cohen, A. (2000), Wavelets in numerical analysis. Handbook of Numerical 
Analysis, vol. VII, P.G. Ciarlet and J.L. Lions, eds.. to appear in 2003 as a 
book "Numerical analysis of wavelet methods", Elsevier, Amsterdam. 

[14] Cohen, A., W. Dahmen and R. DeVore (2000), Adaptive wavelet methods for 
elliptic operator equations — convergence rate. Math. Comp. 70, 27-75. 

[15] Cohen, A., W. Dahmen and R. DeVore (2002), Adaptive wavelet methods for 
operator equations — beyond the elliptic case. Found, of Comp. Math. 2, 203- 
245. 

[16] Cohen, A., W. Dahmen and R. DeVore (2002), Adaptive wavelet methods for 
nonlinear variational problems, preprint IGPM-RWTH Aachen, submitted to 
SIAM J. Num. Anal. 

[17] Cohen, A., S.M. Kaber, S. Mueller and M. Postel (2002), Fully adaptive mul- 
tiresolution finite volume schemes for conservation laws. Math, of Comp. 72, 
183-225. 

[18] Cohen, A. and R. Masson (1999), Wavelet adaptive methods for elliptic 
problems — preconditionning and adaptivity, SIAM J. Sci. Comp. 21, 1006-1026. 

[19] Dahlke, S., W. Dahmen and K. Urban (2001), Adaptive wavelet methods for 
saddle point problems — optimal convergence rates, preprint IGPM- Aachen 

[20] Dahlke, S. and R. DeVore (1997), Besov regularity for elliptic boundary value 



620 



Albert Cohen 



problems, Communications in PDEs 22, 1-16. 
[21] Dahmen, W. (1997), Wavelet and multiscale methods for operator equations, 

Acta Numerica 6, 55-228. 
[22] Dahmen, W., B. GottschHch-Miiller and S. Miiller (2001), Multiresolution 

Schemes for Conservation Laws, Numerische Mathcmatik 88, 399-443. 

[23] Daubechies, I. (1992), Ten lectures on wavelets, SIAM, Philadelphia. 

[24] DeVore, R. (1997), Nonlinear Approximation, Acta Numerica 51-150. 

[25] DeVore, R. and B. Lucicr (1990), High order regularity for conservation laws, 
Indiana J. Math. 39, 413-430. 

[26] Dorfler, W. (1996), A convergent adaptive algorithm for Poisson's equation, 
SIAM J. Num. Anal. 33, 1106-1124 

[27] Eriksson, K., D. Estep, P. Hansbo, and C. Johnson (1995), Introduction to 
adaptive methods for differential equations. Acta Numerica 4, Cambridge Uni- 
versity Press, 105-158. 

[28] Harten, A. (1994), Adaptive multiresolution schemes for shock computations, 
J. Comp. Phys. 115, 319-338. 

[29] Harten, A. (1995), Multiresolution algorithms for the numerical solution of 
hyperbolic conservation laws. Comm. Pure and Appl. Math. 48, 1305-1342. 

[30] Morin, P., R. Nocctto and K. Sicbcrt (2000), Data oscillation and convergence 
of adaptive FEM, SIAM J. Num. Anal. 38, 466 488. 

[31] Pctrushcv, P. (1988), Direct and converse theorems for spline and rational ap- 
proximation and Besov spaces, in Function spaces and applications, M. Cwikel, 
J. Peetre, Y. Saghcr and H. Wallin, eds., Lecture Notes in Math. 1302, Springer 
Verlag, Berlin, 363-377. 

[32] Temlyakov, V. (1998), Best N -term approximation and greedy algorithms. Adv. 
Comp. Math. 8, 249 265. 

[33] Verfiirth, R. (1994), A-posteriori error estimation and adaptive mesh refine- 
ment techniques. Jour. Comp. Appl. Math. 50, 67-83. 



